      subroutine modelo2dcamadas()

      use globais
      implicit none

      double precision::xdonoh,zdonoh
      
      real, dimension(:,:), allocatable:: creal !< Matriz auxiliar, velocidade de propagação em real*4 para impressão em arquivo binário a ser plotado pelas rotinas do Seismic Unix.

      double precision ::xmin,zmin
      integer :: i,j
!velocidade de propagacao no ar = 340m/s
 !velocidade de propagacao na agua = 1450m/s
! massa especifica da agua = 1000

! ! ! ! Modelo da bolinha inspirado no curso do brossier
! !     material 1
!       do j=1,nz
! 	do i=1,nx
! ! 	    c(i,j) = 2250.d0	
! 	    c(i,j) = 4000.d0	
! 	enddo
!       enddo
! 
!       do j=1,nz
! 	zdonoh=   dfloat(j-1)*delta
! 	do i=1,nx
! 	  xdonoh=dfloat(i-Npml)*delta
! 	  if(((xdonoh-Lx0/2.0d0)**2+(zdonoh-Lz0/2.0d0)**2).lt.(200.d0)**2)then
! 	      c(i,j) = 4200.d0	
! 	  endif
! 	enddo
!       enddo
!       cminreal=4000.d0





! ! ! ! Modelo do retangulo inspirado no curso do brossier
! !     material 1
!       do j=1,nz
! 	do i=1,nx
! 	    c(i,j) = 2000.d0	
! 	enddo
!       enddo
! 
!       do j=1,nz
! 	zdonoh=   dfloat(j-1)*delta
! 	do i=1,nx
! 	  xdonoh=dfloat(i-Npml)*delta
! 	  if(xdonoh.gt.150 .and. xdonoh.le.250 .and. zdonoh.gt.100.and. zdonoh.le.200)then
! 	      c(i,j) = 3000.d0	
! 	  endif
! 	enddo
!       enddo
!       cminreal=2000.d0


! ! ! mais um modelo de retangulo 
!     material 1
      do j=1,nz
	do i=1,nx
	    c(i,j) = 2000.d0	
	enddo
      enddo

      do j=1,nz
	zdonoh=   dfloat(j-1)*delta
	do i=1,nx
	  xdonoh=dfloat(i-Npml)*delta
	  if(xdonoh.gt.300 .and. xdonoh.lt.700 .and. zdonoh.gt.300.and. zdonoh.lt.400)then
	      c(i,j) = 4000.d0	
	  endif
	enddo
      enddo
      cminreal=2000.d0

! ! ! ! Modelo inspirado no hu2009, figura 4 do artigo
! !     material 1
!       do j=1,nz
! 	do i=1,nx
! 	    c(i,j) = 2250.d0	
! ! 	    c(i,j) = 1500.d0	
! 	enddo
!       enddo
! 
! !     material 2
!       do j=1,nz
! 	zdonoh=   dfloat(j-1)*delta
! 	do i=1,nx
! 	  if(zdonoh.ge. 100.0) c(i,j) = 1500.d0
! 	enddo
!       enddo
! !     material 3
!       do j=1,nz
! 	zdonoh=   dfloat(j-1)*delta
! 	do i=1,nx
! 	  if(zdonoh.ge. 150.0) c(i,j) = 2000.d0
! 	enddo
!       enddo
! 
! !     material 4
!       do j=1,nz
! 	zdonoh=   dfloat(j-1)*delta
! 	do i=1,nx
! 	  xdonoh=dfloat(i-Npml)*delta
! 	  if(xdonoh.le.120.0 .and. zdonoh.ge.280.0-10.) c(i,j) = 2500.d0
! 	  if(xdonoh.ge.120.0 .and. xdonoh.le.160.0 .and. zdonoh.ge.(-xdonoh+400.0-10.)) c(i,j) = 2500.d0
! 	  if(xdonoh.ge.160.0 .and. xdonoh.le.460.0 .and. zdonoh.ge. 240.0-10.) c(i,j) = 2500.d0
! 	  if(xdonoh.ge.460.0 .and. xdonoh.le.500.0 .and. zdonoh.ge.(xdonoh-220.-10.0)) c(i,j) = 2500.d0
! 	  if(xdonoh.ge.500.0 .and. zdonoh.ge.280.0-10.) c(i,j) = 2500.d0
! 	enddo
!       enddo
! !     material 5
!       do j=1,nz
! 	zdonoh=   dfloat(j-1)*delta
! 	do i=1,nx
! 	  xdonoh=dfloat(i-Npml)*delta
! 	  if(xdonoh.le.240.0 .and. zdonoh.ge.320.0-10.) c(i,j) = 3000.d0
! 	  if(xdonoh.ge.240.0 .and. xdonoh.le.280.0 .and. zdonoh.ge.(xdonoh+80.0-10.)) c(i,j) = 3000.d0
! 	  if(xdonoh.ge.280.0 .and. zdonoh.ge.360.0-10.) c(i,j) = 3000.d0
! 	enddo
!       enddo
!       cminreal=1500.0d0



! ! material 1
!       do j=1,10
! 	  do i=1,Nx
! 	      c(i,j) = 1450.
! 	  enddo
!       enddo
! 
! 
! ! material 2
!       do j=11,int((Nz-Npml+1)*2/3)
! 	  do i=1,Nx/3
! 	      c(i,j) = 2500.
! 	  enddo
!       enddo
!       do j=11,int((Nz-Npml+1)*1/3)
! 	  do i=Nx/3,Nx
! 	      c(i,j) = 2500.
! 	  enddo
!       enddo
! 
! ! material 3
!       do j=int((Nz-Npml+1)*2/3),Nz
! 	  do i=1,Nx/3
! 	      c(i,j) = 5000.
! 	  enddo
!       enddo
!       do j=int((Nz-Npml+1)*1/3),Nz
! 	  do i=Nx/3,Nx
! 	      c(i,j) = 5000.
! 	  enddo
!       enddo
!       
!       cminreal=1450.


! ! material 1
!       do j=1,int((Nz-Npml+1)*2/3)
! 	  do i=1,Nx/3
! 	      c(i,j) = 2000.
! 	  enddo
!       enddo
!       do j=1,int((Nz-Npml+1)*1/3)
! 	  do i=Nx/3,Nx
! 	      c(i,j) = 2000.
! 	  enddo
!       enddo
! 
! ! material 2
!       do j=int((Nz-Npml+1)*2/3),Nz
! 	  do i=1,Nx/3
! 	      c(i,j) = 4000.
! 	  enddo
!       enddo
!       do j=int((Nz-Npml+1)*1/3),Nz
! 	  do i=Nx/3,Nx
! 	      c(i,j) = 4000.
! 	  enddo
!       enddo
!       
!       cminreal=2000.



! ! material 1
!       do j=1,int((Nz-Npml+1)/2)
! 	  do i=1,Nx
! 	      c(i,j) = 1500.
! 	  enddo
!       enddo
! 
! ! material 2
!       do j=int((Nz-Npml+1)/2)+1,Nz
! 	  do i=1,Nx
! 	      c(i,j) = 5000.
! 	  enddo
!       enddo
!       
!       cminreal=1500.



! ! material 1
!       do j=1,31
! 	  do i=1,Nx
! 	      c(i,j) = 2000.
! 	  enddo
!       enddo
! 
! ! material 2
!       do j=32,62
! 	  do i=1,Nx
! 	      c(i,j) = 3000.
! 	  enddo
!       enddo
! 
! ! material 3
!       do j=63,Nz
! 	  do i=1,Nx
! 	      c(i,j) = 5000.
! 	  enddo
!       enddo
!       
!       cminreal=2000.








! ! material 1
!       do j=1,int((Nz-Npml+1)/4)
! 	  do i=1,Nx
! 	      c(i,j) = 1500.
! 	  enddo
!       enddo
! 
! ! material 2
!       do j=int((Nz-Npml+1)/4)+1,int(2*(Nz-Npml+1)/4)
! 	  do i=1,Nx
! 	      c(i,j) = 3000.
! 	  enddo
!       enddo
! 
! ! material 3
!       do j=int(2*(Nz-Npml+1)/4)+1,int(3*(Nz-Npml+1)/4)
! 	  do i=1,Nx
! 	      c(i,j) = 4500.
! 	  enddo
!       enddo
! 
! ! material 4
!       do j=int(3*(Nz-Npml+1)/4)+1,Nz
! 	  do i=1,Nx
! 	      c(i,j) = 6000.
! 	  enddo
!       enddo
! 
! 

! ! velocidade uniforme no modelo
!       do j=1,Nz
! 	  do i=1,Nx
! 	      c(i,j) = 4000.
! 	  enddo
!       enddo
!       cminreal=4000.

! densidade uniforme no modelo
      do j=1,Nz
	  do i=1,Nx
	      rho(i,j) = rhocte
	  enddo
      enddo







! ! ! !       allocate(creal(Nx,Nz))
! ! ! !       open(1,file='veloc.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx*Nz)
! ! ! !       read(1,rec=1)((creal(i,j),j=1,Nz),i=1,Nx)
! ! ! !       close(1)
! ! ! !       do j=1,nz
! ! ! ! 	  do i=1,nx
! ! ! ! 	      c(i,j) = DBLE(creal(i,j)) !< para passar da matriz real*4 para double precision
! ! ! ! 	  enddo
! ! ! !       enddo
! ! ! ! 
! ! ! !       open(1,file='densi.bin',status="unknown",form="unformatted",access="direct",recl=4*Nx*Nz)
! ! ! !       read(1,rec=1)((creal(i,j),j=1,Nz),i=1,Nx)
! ! ! !       close(1)
! ! ! !       do j=1,nz
! ! ! ! 	  do i=1,nx
! ! ! ! 	      rho(i,j) = DBLE(creal(i,j)) !< para passar da matriz real*4 para double precision
! ! ! ! 	  enddo
! ! ! !       enddo
! ! ! !       deallocate(creal)
! ! ! ! 
! ! ! !       cminreal = c(1,1)
! ! ! !       do i=1,nx
! ! ! ! 	  do j=1,nz
! ! ! ! 	      if(c(i,j).lt.cminreal)then
! ! ! ! 		  cminreal=c(i,j)
! ! ! ! 	      endif
! ! ! ! 	  enddo
! ! ! !       enddo








! buoyancy e bulk modulus no grid
      do i=1,Nx
	  do j=1,Nz
	      b(i,j) = 1.0d0/rho(i,j)
	      K(i,j) = rho(i,j)*c(i,j)**2	  
	  enddo
      enddo



      return

      end subroutine modelo2dcamadas